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ABSTRACT 

It is generally accepted that astrophysical sources cannot emit synchrotron radiation above 160 MeV 
in their rest frame. This limit is given by the balance between the accelerating electric force and the 
radiation reaction force acting on the electrons. The discovery of synchrotron gamma-ray flares in 
the Crab Nebula, well above this limit, challenges this classical picture of particle acceleration. To 
overcome this limit, particles must accelerate in a region of high electric field and low magnetic field. 
This is possible only with a non-ideal magnetohydrodynamic process, like magnetic reconnection. We 
present the first numerical evidence of particle acceleration beyond the synchrotron burnofF limit, 
using a set of 2D particle-in-cell simulations of ultra-relativistic pair plasma reconnection. We use 
a new code, Zeltron, that includes self-consistently the radiation reaction force in the equation of 
motion of the particles. We demonstrate that the most energetic particles move back and forth across 
the reconnection layer, following relativistic Speiser orbits. These particles then radiate > 160 MeV 
synchrotron radiation rapidly, within a fraction of a full gyration, after they exit the layer. Our 
analysis shows that the high-energy synchrotron flux is highly variable in time because of the strong 
anisotropy and inhomogeneity of the energetic particles. We discover a robust positive correlation 
between the flux and the cut-off energy of the emitted radiation, mimicking the effect of relativistic 
Doppler amplification. A strong guide field quenches the emission of > 160 MeV synchrotron radiation. 
Our results are consistent with the observed properties of the Crab flares, supporting the reconnection 
scenario. 

Subject headings: Acceleration of particles — Magnetic reconnection — Radiation mechanisms: non- 
thermal — ISM: individual (Crab Nebula) 



1. INTRODUCTION 

The maximum energy reached by a charged particle 
in a given astrophysical object is limited by the size of 
the acceleration region (Hillas 1984). If the relativistic 
Larmor radius of the particle R is of order the system 
size L, the particle escapes and is no longer accelerated. 
The maximum energy is then given by ^max ^ qBL, 
where q is the charge of the particle and B is the typical 
magnetic field strength. Radiative losses within the ac- 
celerator decrease this limit (see, e.g., Aharonian et al. 
2002; Medvedev 2003). The maximum energy is then 
set by the balance between the electric acceleration rate 
and the radiative power lost by the particle. In the 
case of synchrotron cooling, this balance leads to a re- 
markable result: the maximum synchrotron photon en- 
ergy emitted by an electron depends only on the ratio 
of the electric field to magnetic field perpendicular to 
the particle's motion, i.e., e™^^ = {dmc^ /4:aF){E/ B±), 

where ap is the fine structure constant and mc^ is 
the rest mass energy of the electron. Hence, under 
ideal magnetohydrodynamic (MHD) conditions where 
E < B, the energy of synchrotron radiation should 
not exceed the fundamental constant Qmc^/AaF ~ 
160 MeV (Guilbert et al. 1983; de Jager et al. 1996; 
Lyutikov 2010; Uzdensky et al. 2011). An electron with 
energy above the radiation reaction limit would lose most 
of its energy in a fraction of a Larmor gyration. It is then 



impossible to have electrons radiating synchrotron radi- 
ation above 160 MeV with classical models of particle 
acceleration, all based on ideal MHD (e.g., diffuse shock 
acceleration), unless the plasma has a relativistic bulk 
motion with respect to the observer, or the electrons are 
the by-product of energetic particle decay. 

Yet, the gamma-ray space telescopes Agile and Fermi 
discovered several powerful gamma-ray flares from the 
Crab Nebula (Tavani et al. 2011; Abdo et al. 2011; 
Balbo et al. 2011; Striani et al. 2011; Buchler et al. 2012; 
Striani et al. 2013), presumably during which PeV elec- 
trons and positrons emit synchrotron radiation well 
above the 160 MeV limit. The most powerful fiare, 
recorded in April 2011, showed clear evidence for syn- 
chrotron emission up to 375 MeV (Buehler et al. 2012). 
This discovery suggests that an extreme and non- 
conventional particle acceleration mechanism is at work 
somewhere in the nebula, unless the emission is sub- 
stantially Doppler-boosted by a factor > 2. However, 
the typical flow velocity in the nebula (about half the 
speed of light, e.g., Hester et al. 2002) and its orienta- 
tion with respect to the observer give a Doppler factor 
of order unity. The precise location of the flare is un- 
known because the nebula is not resolved in gamma rays, 
but it should not originate very close to the pulsar be- 
cause of lack of pulsations. The < 8 hours flux-doubling 
timescale observed during the flare indicates that a tiny 
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fraction of the nebula is involved. The flaring region ra- 
diates about 30 times more than the quiescent emission 
above 100 MeV, which represents up to 1% of the pulsar 
spin-down power (Buehler et al. 2012). To explain parti- 
cle acceleration and emission within the overall duration 
of the flares, ranging from a few days to a few weeks, 
the magnetic field should be of order a few mG, i.e., 
much more intense than the average ~ 200//G tradition- 
ally inferred from spectral modeling (Horns feAharonian 
2004; Meyer et al. 2010). Simultaneous observations in 
radio, near-infrared, optical. X-rays and in TcV gamma- 
rays, were not able to detect a solid counterpart to the 
flares (see Weisskopf et al. 2013, and references therein), 
suggesting that the emitting particle spectrum is very 
hard, perhaps monocnergetic. It is very difficult to rec- 
oncile these puzzling features with classical models of 
particle acceleration and models of pulsar wind nebu- 
lae (Rees & Gunn 1974; Kennel & Coroniti 1984, and see 
Kirk et al. 2009; Arons 2012 for recent reviews). 

Various models have been proposed to solve the Crab 
flares mystery. Several studies invoke a relativistic 
Doppler boosting of the flaring region by a factor of 
a few. It was proposed that a mildly relativistic flow 
could be achieved close to the pulsar wind termina- 
tion shock (Komissarov & Lyutikov 2011; Lyutikov et al. 
2012; Bednarek & Idee 2011), in relativistic reconnection 
events within the nebula (Clausen-Brown & Lyutikov 
2012), in a magnetized flow at the base of the Crab 
jets (Lyubarsky 2012), or in knots of energetic particles 
(Yuan et al. 2011). The dissipation of the striped pulsar 
wind structure through the shock (Petri & Lyubarsky 
2007; Sironi & Spitkovsky 2011) could also generate 
rapidly fluctuating magnetic field on small scales re- 
sulting in synchrotron gamma-ray flares (Bykov et al. 
2012). If the magnetic turbulence in the nebula oc- 
curs on length scale shorter than the synchrotron pho- 
ton formation length mc^/eB, then the particles could 
radiate jitter radiation (Medvedev 2000), with typi- 
cal energy greater than the classical synchrotron limit 
(Teraki & Takahara 2013). Alternatively, particle accel- 
eration could occur in regions of strong coherent electric 
field, in twisted toroidal fields (Sturrock & Aschwandcn 
2012), or in magnetic reconnection sites within the neb- 
ula (Uzdensky et al. 2011; Cerutti et al. 2012a). 

Magnetic reconnection offers natural locations (within 
the diffusion region where the magnetic field is small 
and reverses) in which the electric field can exceed the 
magnetic field. In principle, it is possible to accelerate 
particles above the classical radiation reaction limit at 
these sites (Kirk 2004). Uzdensky et al. (2011) demon- 
strated that the highest energy particles arc trapped and 
focused towards the reconnection layer midplanc^, fol- 
lowing the relativistic analog of Speiser orbits (Speiser 
1965). Once deep inside the layer, the particles are sub- 
ject to weak radiative losses, but strong coherent elec- 
tric field. The layer acts as a linear accelerator, and 
the maximum energy of the particles is then limited 
just by the total electric potential drop along the layer, 
i.e., Smax eEL, where L is the length of the layer. 
Using relativistic test-particle simulations, we found in 

^ A similar phenomenon occurs in man-made accelerators, where 
gradients of magnetic fields are generated to confine and focus par- 
ticle orbits, see e.g., Courant &; Snyder (1958). 



Cerutti et al. (2012a) that magnetic reconnection could 
generate a quasi-monoenergetic beam of particles above 
the radiation reaction limit in the Crab Nebula. 

In this paper, we re-examine particle acceleration 
in ultra-relativistic pair plasma reconnection, using 2D 
particle-in-cell (PIC) simulations. Such simulations cap- 
ture self-consistently the time-evolution of fields and par- 
ticles at the kinetic level. For this study, we developed a 
new relativistic PIC code, called Zeltron, that includes 
the radiation reaction force on the particles. Our first 
objective is to investigate ab initio whether reconnec- 
tion can accelerate particles above the radiation reaction 
limit, under realistic physical conditions. Our second ob- 
jective is to study the radiative signature of such acceler- 
ation following Cerutti et al. (2012b), in order to explain 
all observed features of the Crab fiares. Cerutti et al. 
(2012b) included optically thin synchrotron radiation as 
a tracer, but this calculation was not self-consistent be- 
cause it did not treat radiation reaction effects on par- 
ticle motion. Section 2 describes the main capabilities 
of Zeltron, with an emphasis on the radiation reaction 
force. Section 3 gives the initial setup of the simulations 
performed in this study, and Section 4 presents the main 
results on particle acceleration and radiation. We discuss 
our findings in the context of the Crab gamma-ray flares 
in Section 5. Section 6 summarizes the main results of 
this work. 

2. THE PARTICLE-IN-CELL CODE ZELTRON 

Zeltron is a new three-dimensional, parallel (domain 
decomposition with MPI), relativistic, electromagnetic 
particlc-in-ccU code developed independently from other 
existing codes (for a review about PIC methods, sec e.g., 
Pritchett 2003; BirdsaU & Langdon 2005). Zeltron fol- 
lows an explicit finite-difference scheme on a Cartesian 
grid, with a time step At set at a fraction of the maxi- 
mum stable step determined by the Courant-Friedrichs- 
Lewy condition, i.e., c/S.t < cAicFL = (1/Ax^ -I- 1/Ay^-f 
1/Az^)~^/^, where c is the speed of light. Ax, Ay, and 
Az are the minimum grid spacing in the x-, y-, and z- 
directions. The code uses the Yee algorithm (Yee 1966) 
to solve the time-dependent Maxwell's equations, given 

by 

— =cVxB-47rJ (1) 



where E is the electric field, B is the magnetic field, and 
J is the current density. This algorithm has second-order 
error in space and time and ensures that V • B = at 
any instant of the simulation (to the computer round- 
off accuracy). The numerical scheme does not, however, 
satisfy the Maxwell-Gauss equation V • E = Airp exactly, 
where p is the charge density. The electric field has to 
be corrected every time step by the small amount JE, 
obtained by solving Poisson's equation (Pritchett 2003), 
i.e., 

(Scj)) = - (47rp - V • E) , (3) 

and (5E = — V(5(/). The Poisson solver implemented in 
the code utilizes an iterative Gauss-Seidel method (with 
5-points in 2D, and 7-points in 3D). 
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The main novelty of Zeltron, compared to most PIC 
codes^, is the ability to take into account the effect of 
the radiation reaction force on the motion of the parti- 
cles. The (non-covariant) equation of motion of a single 
particle is given by the Lorentz-Abraham-Dirac equation 
(Landau & Lifshitz 1975) 



du 

mc- = q[Er 



u X Bj 

7 



(4) 



where u = 7v/c is the four- velocity of the particle, 

7 = (l — /c^) is the associated Lorentz factor, and 
g is the radiation reaction force. The fields at the loca- 
tion of the particle Ej and Bj, are linear interpolations 
of the fields E and B known at the grid nodes. Within 
the framework of classical electrodynamics, the radia- 
tion reaction force is obtained from the Landau-Lifshitz 
equation (Landau & Lifshitz 1975), valid as long as the 
product "fB is much smaller than the quantum critical 
magnetic field i?QED = m^c^ /eh = 4.4 x 10^"^ G, where 
e is the elementary electric charge, and h = hjlix with 
h the Planck constant. In the ultra- relativistic regime 
(7 ^ 1), the radiation reaction force can simply be ex- 
pressed as a continuous friction force opposite to the par- 
ticle's direction of motion (see discussion in Cerutti et al. 
2012a), i.e.. 



' 2 



Ei 



u X Bj 

7 



u Ei 

7 



(5) 



where — e^/mc^ is the classical radius of the elec- 
tron. Following Tamburini et al. (2010), Zeltron uses a 
modified Boris algorithm (Birdsall & Langdon 2005) to 
solve the equation of motion with the radiation reaction 
force given in Eq. (5). For an alternative implementa- 
tion of the radiation reaction force in PIC codes, see also 
Sokolov et al. (2009); Capdessus et al. (2012). 

The code uses linear interpolation to deposit the 
charges and currents generated by each particle at the 
nodes of the computational grid, and computes the 
charge and current densities for Maxwell's equations. 
The code assigns variable weights to the macro-particles 
to model particle density gradient. Zeltron does not 
strictly conserve the total energy. For the purpose of 
this study, we ran Zeltron successfully on thousands of 
cores on the Kraken supercomputer^ and on the Univer- 
sity of Colorado Janus and Verus supercomputers. We 
find excellent agreement between Zeltron and the well- 
tested PIC code Vorpal (Meter & Gary 2004), in the 
limit where radiative losses are neglected (the radiation 
reaction force is currently not implemented in Vorpal). 

3. SETUP OF THE SIMULATIONS 

^ Several PIC codes used in the laser-plasma interaction 
community (see e.g., Zhidkov et al. 2002; Sokolov ct al. 2009; 
Tamburini et al. 2010; Capdessus et al. 2012) do include the ra- 
diation reaction force, in preparation for future radiation pressure- 
dominated plasma experiments with the new generation of ultra- 
intense lasers. We note also that Jaroschok & Hoshino (2009) 
present the first study of radiation-dominated rcconncction in the 
relativistic regime, using PIC simulations with radiation reaction 
force. 

^ National Institute for Computational Sciences 
(www.nics.tennessee.edu/). 



We present a series of simulations of ultra-relativistic 
electron-positron pair plasma reconnection with 
radiation reaction force in 2.5D (the fields de- 
pends on 2 coordinates, but the motion of parti- 
cles is 3D), using Zeltron. The initial setup is 
very similar to our previous study (Cerutti et al. 
2012b), which is standard in such reconnection 
simulations (Zenitani & Hoshino 2001, 2007, 2008; 
Jaroschek et al. 2004b; Bessho & Bhattacharjee 2007, 
2012; Daughton & Karimabadi 2007; Petri & Lyubarsky 
2007; Yin et al. 2008; Jaroschek & Hoshino 2009; 
Liu et al. 2011; Sironi & Spitkovsky 2011; Kagan et al. 
2012). The computational domain is a rectangle of size 
Lx X iy, where Lx.y is the length of the box in each 
direction. It contains two anti-parallel relativistic, flat, 
Harris current sheets (Kirk & Skjasraasen 2003) in the 
a;2;-plane, which enables us to set periodic boundary 
conditions in the x- and y- directions. The reconnecting 
magnetic field, i3x, reverses across each layer along the 
y-direction, and is given by 



—Bo tanh 
Bn tanh 



S 



if y < Ly/2 
if y > Ly/2 



(6) 



where Bq is the initial upstream reconnecting field 
strength, and S is the initial layer half-thickness. The 
strength of the uniform guide field component (i?z) is a 
free parameter that varies from to Bq in our simula- 
tions (see Section 4.5). There is no electric field initially, 
E = 0. We apply an initial perturbation to the mag- 
netic field (10% maximum amplitude in the magnetic 
flux function, see Fig. 1), in order to force the onset of 
reconnection at the beginning of the simulation, thereby 
reducing the computing time. If there is no perturba- 
tion, we flnd that reconnection eventually happens spon- 
taneously in the simulation because of the high numerical 
noise inherent to PIC codes. 

A population of relativistic thermal pairs, following 
a Maxwell- Jiittner distribution and concentrated in the 
current layers, balances the upstream magnetic pressure 
and carries the initial current in the ±2;-directions given 
by the bulk motion of these particles, drifting at a ve- 
locity Wdrift/c = /3drift — 0.6. The density profile of the 
drifting particles is 



no 



no 



cosh 



cosh 



(^) 

j/-3£y/4 



if y < Ly/2 
-2 , (7) 

if y > Ly/2 



with 



no 



fcTdHft(l-/3^Hft) 



1/2 



47re2/3, 



(8) 



drift 



where k is the Boltzmann constant, and Tdrift is the 
temperature of the drifting particles defined in their co- 
moving frame. In the simulation, the drifting particles 
are injected uniformly throughout the box, but with vari- 
able weight according to Eq. (7). The purpose of this pro- 
cedure is to represent low density distributions with less 
noise. In addition, we fill the box with a uniform density 
ribg = O.lno of isotropic non-thermal ultra-relativistic 
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particles (see explanation below), such that 

^^bg ^ f Kj-P if 71 < 7 < 72 (g-. 

1 otherwise ' ^ ' 

p is the index of the power-law, and X is a normalization 
constant. 

We initialize the simulations with physical parameters 
consistent with the conditions thought to exist in the 
flaring region of the Crab Nebula, although our results 
are general and scalable. The energy and spatial scales 
of the problem are obtained from the reconnecting field 
strength Bq, of order a few milliGauss during Crab flares. 
Following Cerutti et al. (2012a), we set Bq = 5 mC. The 
fiducial radiation-reaction-limited energy of the particles, 
given by the balance of the electric force with the radia- 
tion reaction force for E = Bq, is then equal to 

^'^^ = (2^1^) " , (10) 

where i?5mG = Bq/5 mG, allowing electrons with ener- 
gies up to in the PeV range. The Larmor radius asso- 
ciated with such a particle is i?rad = Ira.d'mc^ / cBq k, 
4.6 X 10^''^ cm, or 4.25 light-hours. The maximum energy 
expected to be reached by the particles in the simula- 
tion (in the absence of radiative losses) is limited by the 
electric potential drop along the layer, which is propor- 
tional to the system size L and the dimensionless recon- 
ncction rate /3roc < 1, such that 7niax ^ P-ccc^^BqL / mc^ . 
To observe particle acceleration above the radiation re- 
action limit (Eq. 10), the system size L must be bigger 
than i?rad- -Rrad Hiust also be large compared with the 
smallest spatial scales in the system, given by the Lar- 
mor radius of the lowest energy particles injected ini- 
tially. We set the minimum Lorentz factor of the back- 
ground particles to 71 = 4 x lO*" with a power-law index 
of j3 = 2 extending up to 72 = 4 x 10^ < 7rad, and an 
ultra-relativistic temperature /cTdrift/wc^ = 4 x 10'' for 
the drifting particles. In reality, the typical energy of 
the background particles in the Crab Nebula is proba- 
bly much lower, of order 71 ^ 10^-10^. However, such 
a broad range in energy cannot be reached with the 
current computing power available, even for 2D simu- 
lations. Our simulations model only the high-energy end 
of the Crab Nebula spectrum. The initial —2 power- 
law of the electrons approximately represents the purely 
non-thermal background particles emitting the quiescent 
emission of the Crab Nebula, inferred from spectral mod- 
eling (Horns feAharonian 2004; Meyer et al. 2010). This 
assumes pre-acceleration of the particles in the nebula, 
possibly by shock acceleration or by other reconnection 
events. 

In the following, we express all spatial scales in terms of 
the initial minimum Larmor radius pi = jimc^ / cBq « 
1.4 x 10^^ cm, and timescales with respect to the cor- 
responding inverse Larmor frequency oj^^ = pi/c w 
455 seconds. We perform all the simulations with = 
Ly = 500pi, which corresponds to 7x 10^'' cm or 2.7 light- 
days. For a typical value of /3,.oc = 0.1-0.2, we have 
7max £ 5O-IOO71 « 2-4 X 10^, allowing us to see particles 
with gamma above 7rad ~ 1.3 x 10^. The other relevant 
quantities of the plasmas are the electron skin-depth in 
the layer de — (fcT(irift/47rnoe^)^/^ w l.Spi, the initial 



layer thickness 6 = 2fcTdrift(l - PlittV^V l^d.meBo ~ 
2.7 pi, and the upstream magnetization parameter a ~ 
i?Q/47rnbg7i mc^ ss 16. The computational grid is Carte- 
sian and uniform, composed of 1440^ cells, giving a spa- 
tial resolution of w 3 cells per pi or w 8 cells per 5. We 
note that the simulation becomes unstable at late times 
if At = AtcFL, only when radiative cooling is strong 
and if there is no guide field. In particular, the elec- 
tric field oscillates between two time steps leading to an 
artificial increase of radiative losses. We find that set- 
ting the time step to At = 0.3AtcFL is good enough 
to quench the development of this numerical instability. 
The time resolution is then At ~ O.QToJi^. The total 
initial number of particles per cell is 100 (25 electrons 
and 25 positrons for each population). Table 1 gives all 
the physical and numerical parameters and their values 
chosen for this study. 

TABLE 1 

Physical and numerical parameters set up in all the 
simulations reported in this work. 



Physical parameters 


Set values 


Bo 


5 mG 


7rad 


1.3 X 10*' 


Pi 


1.4 X 10^^ cm 


dc/pi 


1.8 




2.7 




455 s 




4 X 10^ 


/3drift 


0.6 


"bg/no 


0.1 


a 


16 


71 


4 X 10^ 


72 


4 X 10* 


P 


2 


Numerical parameters 


Set values 


pi /Ax 


3 




3 


At X u)\ 


0.07 


Particles / cell 


100 


Lyi/pi 


500 


Ly/Pl 


500 



4. RESULTS OF THE SIMULATIONS 

This section presents the results based on a set of simu- 
lations of size {L^xLy) = (500/9i)^ to investigate the gen- 
eral properties of reconnection under strong synchrotron 
cooling conditions and particle acceleration above the 
radiation reaction limit (Section 4.1-4.3); we compare 
with an identical simulation without the radiation re- 
action force. In Section 4.4, we investigate the variabil- 
ity patterns (spectra, light curves, power-spectra) of the 
high-energy radiation emitted by the layer, and we show 
that these features are robust and reproducible, using 10 
identical simulations with the radiation reaction force. 
The statistical variations in this sample of simulations 
originate only from the initial random positions and ve- 
locities of the particles in the box. We present also a set 
of 4 simulations to study the effect of the guide field on 
particle acceleration above the radiation reaction limit 
(Section 4.5). 

4.1. Time evolution of reconnection 
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Fig. 1. — Plasma density (color-coded) normalized to the initial drifting particle density no, and magnetic flux function iso-contours 
tracing magnetic field lines (white solid lines), shown at tuji = (top-left), tuji = 198 (top-right), tuji = 309 (bottom-left), and toji = 662 
(bottom-right). 



At the beginning of the simulation, the system remains 
quasi-static"*^ for about tuji < 110, before the layers be- 
come unstable to tearing modes and break up into several 
islands of closed magnetic field loops filled with plasma 
("magnetic islands" or "plasmoids", see Fig. 1). This in- 
stability produces multiple secondary X-points between 
the plasmoids, where the magnetic field reconnects and 
where the electric field, mostly along the ±2;-directions 
{E fa Ez) is intense, of order Bq, leading to strong par- 
ticle acceleration (see Fig. 2, bottom-left panel). In ad- 
dition, the formation of islands induces strong bipolar 
reconnected magnetic field By ^ Bq concentrated at the 
edges of islands between two X-points (see Fig. 2, middle- 
right panel) . The magnetic tension of this field drives the 
rcconncction outfiow by pushing the plasma away from 
X-points into the direction of magnetic islands. Later on, 
magnetic islands become rounder and merge with each 
other to form bigger islands, until there is only a single 
big island and a single X-point left (per layer) at the end 
of reconnection (at ttoi > 450, Fig. 1; the end of the 
simulation does not reach the full saturated state). This 
peculiar symmetric final state arises from the choice of 
double-periodic boundary conditions. 

* This phase is significantly shortened by the initial perturbation 
applied to the magnetic field. 



We find that the rcconncction rate increases by a fac- 
tor of 2 when the radiation reaction is self-consistently 
included, /3ioc ~ 0.3. Because the initial pressure balance 
is approximately maintained during the simulation, the 
layer compresses due to high synchrotron energy losses 
inside the layer, leading to an increase in the reconnection 
rate (Uzdensky & McKinney 2011). The magnetic is- 
lands tend to be smaller and denser with radiative losses. 
Apart from this, we do not find any qualitative difference 
in the overall time-evolution of the reconnection process, 
with or without the radiation reaction force. Fig. 3 shows 
the time evolution of the different energy components in 
the system, i.e., the energy of the fields, of the particles 
and of the radiation. About 63% of the initial magnetic 
energy is dissipated and converted to energetic particles 
and radiation. The total energy is conserved to within 
about 5% throughout the simulation. This moderate er- 
ror comes from the overestimation of the synchrotron 
energy losses due to the high-frequency fiuctuating elec- 
tric field described in Section 3. If the radiation reaction 
force is neglected, or if there is a non-zero guide field, the 
total energy is very well conserved with less than 0.2% 
error. 

To summarize, the overall time evolution of reconnec- 
tion can be divided schematically into three main phases 
shown in Fig. 3: 
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Fig. 2. — Electric (loft) and magnetic (right) field components normalized to Bq. The x- (i?x, Bx, top panels), y- {Ey, By, middle panels), 

198, zoomed in on the lower layer islands. White solid lines represents 



and z- (_Ez, Bz, bottom panels) components arc shown at tui 
magnetic field lines. 
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Fig. 3. — Time evolution of the total magnetic ("i?mag", blue 
solid line) and electric field energies ("i?oicc"i multiplied by a fac- 
tor 20 for readability, red dot-dashed line), the kinetic energy of 
all the particles ("-Ekin"i green dashed line), and the total en- 
ergy lost through synchrotron cooling by the particles ("-Ei.ad"j 
purple dotted line). The three main phases of reconnection (1. 
"Quasi-static", 2. "Plasmoid-dominated" , and 3. "Saturated") 
are roughly delimited by vertical dotted lines. 



Phase 1. Quasi-static state where the electric field 
builds up smoothly. 

Phase 2. Plasmoid-dominated reconnection. This 
is the most active and bursty period of magnetic 
dissipation, and there is a strong competition be- 
tween particle acceleration and cooling. 



• Phase 3. Saturated state, where particle cooling 
dominates. The particles and the magnetic field 
are approximately in cquipartition. 

4.2. Particle and synchrotron radiation spectra and 
anisotropics 

Fig. 4 presents the energy distributions of all the back- 
ground particles (7^dN/c?7, top panel) and their instan- 
taneous, transparent synchrotron radiation {vFu, bottom 
panel) at tuji = 0, 110, 220, 397, and 662. The con- 
tribution from the drifting particles is not shown here 
because of their small number compared with the back- 
ground particles. We assume that the photons are emit- 
ted continuously (valid if <^ Sqed) and tangentially 
to the particle's orbit (valid if 7 ^ 1). Photons do 
not interact with the plasma (optically thin approxima- 
tion), hence there is no need to solve the full radiative 
transfer equation. Fig. 4 shows the emergence of a high- 
energy tail of particles, whose maximum energy increases 
with time during the most active phase of reconnection, 
tuji < 450 (i.e., during the plasmoid-dominated phase). 
During the early stage of reconnection, the initial —2 
power-law extends to higher energies, and cuts off ex- 
ponentially. The spectrum marginally extends at lower 
energies. At tuii > 220, there are some particles that 
are accelerated above the radiation-reaction limit energy 
7rad « 1.3 X 10^. We find that these extremely ener- 
getic particles are accelerated at X-points, surfing mul- 
tiple times across the reconnection layer, following rela- 
tivistic Speiser orbits. The details of the Speiser acceler- 
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Fig. 4. — Particle {■y'^ dN / d'y , top) and photon spectral (f-Fi/, bot- 
tom) energy distributions averaged over all directions, at tuii = 
(dotted line), 110, 220, 397, and 662 (dashed line). The red verti- 
cal dot-dashed line in the top panel shows the classical radiation- 
reaction-limited energy 7i.ad defined in Eq. (10), and the corre- 
sponding 160 MeV synchrotron photon energy limit in the bottom 
panel. The thick black 3-dot-dashed line marked "no rad." in 
both panels shoves the final quasi-steady energy distributions at 
tuii = 662, if there arc no radiative losses in the simulation. 

ation mechanism are described below in Section 4.3. At 
tui = 397, we find that about 0.03% of all the particles 
are above 7rad, and represent about 0.39% of the total 
kinetic energy. These particles are responsible for the ex- 
cess of synchrotron radiation above 160 MeV (see Fig. 4), 
which represents about 4.3% of the total isotropic radia- 
tive power. After tuji = 450, the high-end of the spec- 
trum contracts to lower energies, and very few particles 
above 7i-ad survive at tuji > 662. We attribute the disap- 
pearance of the most energetic particles to synchrotron 
cooling. In the final saturated state, most particles are 
located within the big islands where they cool progres- 
sively. There is little acceleration and the magnetic field 
remains strong, of order Bq- If the radiation reaction 
force is neglected, the high-energy component of the 
spectrum does not evolve once established (for tuii > 550, 
see Fig. 4), and extends up to 7max ~ 3 x 10^. An iden- 
tical simulation performed without radiative losses with 
Vorpal gives very similar results. 

We investigate the angular distribution of the parti- 
cles' velocities, as a function of their energy. In agree- 
ment with our previous study (Cerutti et al. 2012b), we 
find a pronounced energy-dependent anisotropy of the 
particles and their synchrotron emission, increasing with 
energy. Fig. 5 illustrates the strong anisotropy of the ex- 
pected synchrotron radiation, as a function of the pho- 
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Fig. 5. — Encrgy-resolvcd angular distribution of the synchrotron 
radiation fiux d{uF^)/dQ/dei emitted at tuJi = 397, using the 
Aitoff projection. Each panel shows the angular distribution of ra- 
diation in a different photon energy band: 1 MeV < ei < 1.2 MeV 
(top), 12.6 MeV < ei < 14.5 MeV (middle), and 155.7 MeV < 
ei < 179.0 MeV (bottom). Fluxes are normalized to the maximum 
value in each band. The solid angle covered by half of the fiux and 
normalized by 47r, f25o/47r, is given below each panel. The black 
square box indicates the direction where the anisotropic spectra 
arc shown in Fig. 6. 

ton energy. Following Cerutti et al. (2012b), we use the 
spherical angles (f> (latitude) and A (longitude) to study 
the angular distributions. The latitude varies between 
—90° and +90°, and the longitude varies between —180° 
and -1-180°. A radial unit vector has the coordinates 
X = cos sin A, y — sin cj), z = cos(/)COsA. At = 397, 
we find that half of the > 160 MeV radiative flux is con- 
centrated into less than 4% of the total solid angle Atk. 
The high-energy beam of radiation is concentrated in the 
mid-plane (xz-plane, = 0°), preferentially towards the 
zbx-directions ((/> = 0°, A = ±90°). This result can be 
explained by the deflection of the particles' trajectories 
from the ±z-directions (</) = 0°, A = 0°, ± 180°, along 
which the particles are accelerated by E^) to the ±a;- 
directions by the reconnected fleld. The direction of the 
beam is changing with time, wiggling around the plane of 
the reconnection layer during the active phases of recon- 
nection (tujx < 450). The beam broadens and stabilizes 
along the z-direction at later times. 

The strong anisotropy of the emitted radiation leads 
to an apparent boosting of the flux seen by an observer 
looking in the direction of the beam (the so-called "ki- 
netic beaming", see Cerutti et al. 2012b). The energy 
distribution of the particles pointing in the direction 
A = -f-70°, = 0° (indicated by the black box in Fig. 5) 
within the solid angle AJ7/47r sa 3 x 10"'^ is substantially 
harder than the isotropic one, piling up at 7niax ~ 7rad 
(Fig. 6, top panel). In this case, the particles with 
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Fig. 6. — Particle {■y'^dN/df, top) and spectral {vF,^, bottom) 
energy distributions at tui\ = 397. The solid lines labeled "ANIS." 
give the distributions as seen by an external observer looking at the 
direction A = +70°, <^ = 0° (i.e., in the plane of the reconnection 
layer, close to the +a;-direction), within a solid angle AQ/Att Ri 
3 X 10"'^ (shown by the blacl< box in Fig. 5). The dashed lines are 
the isotropically averaged distributions labeled "ISO." , as shown in 
Fig. 4 for comparison. The red vertical dot-dashed line marks the 
classical radiation-reaction-limited energy 7i.ad in the top panel, 
and the corresponding 160 MeV synchrotron photon energy limit 
in the bottom panel. 

Lorentz factors above 7rad account for about 0.5% of 
the particles and 8.5% of the energy. Similarly, the ap- 
parent spectral energy distribution is dominated by the 
high-energy radiation, peaking at around 100 MeV with 
about 20% of the radiative power above 160 MeV. The 
spatial distribution of high-energy particles is strongly 
inhomogeneous (see also Cerutti et al. 2012b). Fig. 7 
shows the spatial distribution of a sample of high-energy 
particles with 7 > 5 x 10*, at tcui = 309. The ener- 
getic particles are clustered into compact bunches within 
the reconnection layer and magnetic islands. The pro- 
nounced inhomogeneity and anisotropy of the extremely 
energetic particles above 7rad, and the associated radia- 
tion above 160 MeV, are key elements in explaining the 
Crab gamma-ray flares (see Section 5). 

4.3. Relativistic Speiser orbits 

In this section, we examine in detail the acceleration 
mechanism of the most energetic particles with 7 > 7rad • 
We follow the trajectories of 20, 000 particles, picked 
randomly and uniformly throughout the box at t = 0. 
In this sample, there are about a dozen particles reach- 
ing a maximum energy above 7rad- Fig. 8 shows three 
typical particle trajectories projected in the yz-plane 
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Fig. 7. — Spatial distribution of a random sample of 234 high- 
energy particles with 7 > 5 X 10®. The red squares show the 
locations of each these particles in the xjz-plane at tuii = 309. The 
solid lines are the magnetic field lines. 

and in the xy-plane, as well as the time history along 
the particle trajectory of the particle Lorentz factor 7, 
the relative strength of the electric force and the radia- 
tion reaction force, and the synchrotron critical energy 
Csync = iheB / Airmc. 

These extremely energetic particles systematically fol- 
low the same simple time history, which can be decom- 
posed into four main phases: 

• 1. Drifting: The particle is initially located up- 
stream and is well magnetized where ideal MHD 
holds {E < B±). The particle moves together with 
the field lines towards the layer, gains little energy 
and does not radiate much. This phase ends at the 
reference point "in" in Fig. 8 when the particle gets 
inside the layer. 

• 2. Linear acceleration: This phase starts when 
the particle reaches the reconnection layer, where 
it is no longer magnetized. The strong electric field 

accelerates the particle almost linearly along 
the 2-direction. In addition, the magnetic field 
reversing across the layer confines the particle to- 
wards the layer mid-plane. The particle follows a 
relativistic Speiser orbit (Speiser 1965; Kirk 2004; 
Uzdcnsky et al. 2011), whose meandering width j/m 
(maximum distance of the particle from the neutral 
sheet) decreases with time as the particle energy 
increases. The particle is then confined closer and 
closer to the layer mid-plane, where the perpendic- 
ular magnetic field strength and hence the radia- 
tion reaction force decrease, while the electric force 
keeps on accelerating the particle. During this pro- 
cess, the particle's trajectory is also significantly 
deflected towards the zbx-directions by the recon- 
nected field By. This phase begins at the reference 
point "in" and ends when the particle reaches its 
maximum energy, marked by the red cross in Fig. 8. 
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Fig. 8. — Top panels: Three typical trajectories (projected onto the yz- and xy-planes: note the compressed y-axis) of high-energy 
particles accelerated above the radiation reaction limit, 7radi inside the rcconncction layer. The gray bands indicate the initial thickness 
of the layer. The green triangle marks the position of the particle at t = 0, while the red diamond marks the position of the particle at 
the end of the simulation at tLO\ = 662. Bottom panels: Time evolution along the particle trajectory of (from top to bottom) the Lorentz 
factor, the strength of the electric force Fc(t) (purple solid line) and of the radiation reaction force g{t) (blue dotted line), and the critical 
synchrotron photon energy esync(t)- The black squares indicate two reference points along the particle orbit where the particle moves "in" 
and "out" of the layer. The crosses mark the maximum values of 7 and esync reached by the particle; these values are also printed nearby. 



• 3. Ejection and emission: While it remains 
deep inside the layer, the particle is accelerated 
above the radiation reaction limit 7rad- The parti- 
cle is then ejected from the layer when it encoun- 
ters a magnetic island (in most cases the final big 
magnetic island towards the end of the simulation) 
and feels a sharp increase of the radiation reaction 
force. The particle loses most of its energy in a 
fraction of a Larmor cycle and emits synchrotron 
photons above 160 MeV. This phase happens just 
before and after the particle exits the layer (around 
the reference point "out" in Fig. 8). 



• 4. Cooling: The particle is back in a region 
where ideal MHD conditions apply and cools pro- 
gressively. The particle does not experience any 
significant acceleration once the saturated state of 
rcconncction is reached. 

This simple picture highlights the distinction between 
the acceleration zone (inside the layer) and the > 
160 MeV synchrotron radiating zone (upstream, or in- 
side magnetic islands). The acceleration zone is of order 
the system size ^ L^, while the radiating zone is a frac- 
tion of the Larmor radius ^ ^^^^^mc^ / eB^. The acceler- 
ation and focusing mechanisms described above in phase 
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Fig. 9. — Evolution of the particle's mid-plane crossing angle 
6o with the particle's Lorcntz factor 7, for a representative sam- 
ple of 8 high-energy particles accelerated via the Speiser mech- 
anism. The green triangle and the red diamond mark the first 
and the last crossing of the particle through the layer mid-plane. 
The particles shown here undergo between 4 and 9 crossings before 
they are kicked outside the layer. The arrow along each particle's 
path indicates the direction of increasing time. The power-laws 
of index —2/3 (dashed lines) and —3/2 (dotted lines) are analyti- 
cal solutions of relativistic Speiser orbits found by Uzdensky et al. 
(2011). The vortical dot-dashed line shows the radiation reaction 
limit Lorcntz factor 7i-ad (Eq. 10). 

2 agree surprisingly well with our previous test-particle 
simulations (Uzdensky et al. 2011; Cerutti et al. 2012a), 
despite the simplistic assumptions on the fields used in 
these studies (prescribed and static). Uzdensky et al. 
(2011) predicted a relationship between 7 and the an- 
gle 9q between the particle's velocity vector and the 
layer mid-plane defined at each crossing, in two extreme 
regimes. If the particle's meandering width is much 
greater than the layer thickness (5, and if radiative losses 
are negligible (i.e., during the first Speiser cycles), then 
16*01 oc 7~^/'^. In contrast, if the particle is deep inside the 
layer and reaches the local radiation reaction limit energy 
'Yrad (defined with the perpendicular field at the location 
of the particle B±^ < Bq, so that 7^^^^ > 7rad) within 
each cycle, then \6o\ oc 7^^/^. Fig. 9 shows the tracks 
followed by a representative sample of 8 high-energy par- 
ticles in the Oq-^ plane (which are not necessarily acceler- 
ated above 7rad)- The mid-plane crossing angle is given 
by ^0 = 7''/2 — arccos(wy/-y/v • v), where v is the three- 
velocity vector of the particle. The agreement with the 
analytical expectations is very good: the particles remain 
between these two power-laws, tending to a —2/3 index 
at low energies and to a —3/2 index at the highest en- 
ergies. This is a robust and clean feature of the most 
energetic particles accelerated and focused through the 
Speiser mechanism. 

4.4. Variability pattern of the >100 MeV emission 

In this section, we investigate the time-dependent ra- 
diation escaping in the -|-a;-direction where most of the 
high-energy radiation is expected (Fig. 5). Fig. 10 
presents the expected synchrotron fiux integrated above 
100 MeV as a function of time, taking into account the 
time delay due to the light crossing time through the 
box. In the case of radiation into the -|-a;-direction, the 
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Fig. 10. — Normalized synchrotron flux emitted by the positrons 
as function of time (given in days, bottom axis, and in light crossing 
time of the system, ci/Lx, top axis) in three photon energy bands: 
1 MeV< ei < 10 MeV (green dashed line), 10 MeV< ei < 100 MeV 
(blue dotted line), and ei > 100 MeV (red solid line). The ra- 
diation received by the observer is going along the -l-x-direction 
(0 = 0°, A = -1-90°) throughout the simulation within a solid an- 
gle Ar2 0.03 Sr. The radiation comes from the bottom layer only. 
The vertical dotted lines delimit the 12 time periods of equal dura- 
tion, used to study spectral variability above 100 MeV in Fig. 12. 

propagation time is given by ipropag = (^x — Xc)/c, where 
Xc is the location of the emitting electron/positron. In 
agreement with Cerutti ct al. (2012b), the high-energy 
radiation is highly variable on timescales much shorter 
than the light crossing time of the layer (< O.lLx/c, or 
< 6 hours). The light curve is composed of multiple 
intense spikes that are nearly symmetric in time. This 
result is a direct consequence of the strong focusing of 
the energetic particles accelerated through the Speiser 
mechanism. The beam of energetic particles is wiggling 
around the reconnection layer and crosses the line of sight 
several times. The bunching of the high-energy particles 
into compact blobs within the layer and within the mag- 
netic islands also contributes to the multiple, powerful 
sub-flarcs in the light curve (Cerutti et al. 2012b). This 
dramatic variability disappears if, instead of considering 
one particular direction, the emission is averaged over all 
directions. Fig. 10 also shows the energy dependence of 
the light curve. The amplitude of the spikes increases 
with the energy of the radiation considered, because 
of the increasing emission anisotropy (Fig. 5). Fig. 11 
presents the resulting power-density-spectrum (PDS) of 
the light curve (given by the squared modulus of the 
Fast-Fourier- Transform), in the three energy bands de- 
fined in Fig. 10. The observed PDS above 100 MeV is 
well-fit by a hard power-law of index —0.50 ± 0.05. At 
lower energies, the best-fit indexes arc — 1.04± 0.13 in 
the 10 MeV< ei < 100 MeV band, and -1.19 ± 0.08 in 
the 1 MeV< ei < 10 MeV band. As expected, the PDS 
slope hardens with increasing photon energy, indicating 
that the highest energy radiation is also the most rapidly 
variable. 

The received spectrum is also highly time variable. 
We decompose the light curve into 12 blocks of duration 
12 hours each (see Fig. 10). Within each period of time, 
we compute the time-averaged synchrotron spectrum re- 
ceived by the observer. We fit the high-energy compo- 
nent above 100 MeV only with a power-law times an ex- 
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Fig. 11. — Powcr-density-spcctrum of the observed light curve in 
the three photon energy bands defined in Fig. 10; 1 MeV< ei < 
10 MeV (green fines, top), 10 MeV< ei < 100 MoV (blue Unes, 
middle), and ei > 100 MeV (red lines, bottom). Tfie solid lines 
sfiow the PDS averaged over the statistical sample of 10 simula- 
tions. The dotted lines give the PDS of the light curve shown in 
Fig. 10 only. The frequency / = 1/t is in day~^ (and in units of 
c/Lx, top axis), ranging from the inverse of the total duration of 
the light curve, i.e. fa 1/6 day~^, to the inverse of the time res- 
olution of the light curve At],. ^ 30 day~^. The red dashed lines 
are best-fit power-laws of the mean power-density-spectra, with 
indexes shown above each line. 

ponential cut-off, vF^, = i^^^e" exp (— ei/ccut), where the 
free parameters are the normalization constant Kf:, the 
spectral index a and the cut-off energy Ccut- This analy- 
sis reveals a strong correlation between the cut-off energy 
and the total flux above 100 MeV (see Fig. 12). A power- 
law fit gives i'Fiy{€i > 100 MeV)= /ioo°McV ^'^^^i ^ 
^cift*^"^- To test the robustness of this correlation, we 
performed a series of 10 identical simulations, where only 
the initial positions and velocities of particles vary. We 
find a strong and positive flux/cut-off energy correla- 
tion in all the simulations, using 50 time bins to refine 
our analysis. The power-law index varies from about 
corr « +2 to corr « -1-4, with a mean index corr « -1-3 
(see Tab. 2). We note that the slope depends on the low- 
energy cut-off (100 MeV). We are also able to confirm the 
reproducibility of the power-density-spectrum pattern 
from one simulation to another (Fig. 11). The average 
PDS are best fitted by power laws of indexes — 1.24±0.04 
in the 1 MeV< ei < 10 MeV band, -1.06 ± 0.03 in the 
10 MeV< ei < 100 MeV band, and -0.77 ± 0.03 above 
100 MeV. The uncertainties given here do not represent 
the statistical uncertainties, but the error on the fit only. 

4.5. Effect of the guide field 

All the simulations presented above assumed a zero 
magnetic guide field component = 0. In the general 
case, however, reconnection is not purely anti-parallel 
and there is always a finite guide field. In this section, 
we examine the effect of a moderate and uniform guide 
magnetic field on particle acceleration and radiation, by 
performing a set of 4 simulations with B^/Bq = 0.1, 0.25, 
0.5, and 1, with the physical parameters given in Tab. 1. 
The general time evolution of the reconnection dynamics 
described in Section 4.1 is still valid with a guide field. 
However, we find that the morphology of the particle 



TABLE 2 

Correlation between the energy flux and the high-energy 

PHOTON spectral CUT-OFF ENERGY, FOR A STATISTICAL SAMPLE 
OF 10 IDENTICAL SIMULATIONS WHERE ONLY THE INITIAL POSITIONS 
AND VELOCITIES OF PARTICLES VARY. THE CORRELATION 
COEFFICIENT IS THE BEST-FIT POWER-LAW INDEX AS SHOWN IN 

Fig. 12, I.E., uF,y{€i > 100 MeV)oc ej:"^, using 50 bins in time. 



Simulation 


Best-fit index corr 


1 


3.14 ±0.45 


2 


2.57 ±0.44 


3 


3.26 ±0.41 


4 


2.90 ±0.50 


5 


1.96 ±0.38 


6 


4.07 ±0.43 


7 


3.33 ±0.52 


8 


2.89 ±0.39 


9 


3.05 ±0.37 


10 


3.44 ±0.56 
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e„„i [MeV] 



Fig. 12. — Energy flux of the synchrotron radiation above 
100 MeV {uF^iei > 100 MeV)) as function of the spec- 
tral cut-off energy, ecut, found using the analytical fit uF,y = 
Kytf exp (— ei/ecut)- Each point corresponds to the observed spec- 
trum averaged over the time periods from 1 to 12 defined in the 
light curve in Fig. 10. The periods "1" and "2" do not appear in this 
plot because there is no high-energy flux at these times. The arrow 
shows the path followed over time by the synchrotron spectrum. 
There is a clear correlation between the energy flux and the cut- 
off energy. A power-law fit gives uF^{e-i > 100 MeV)oc e^ut*'*'" '', 
over-plotted here with the red dashed line. The dotted vertical line 
marks ecut = 160 MeV. 

distribution inside the plasmoids is more complex. In 
the presence of a finite guide field, two streams of par- 
ticles form symmetrically on both sides of the layer, in 
the opposite direction to the island motion (Fig. 13, left 
panels). Each arm is produced by the deflection of the 
particles by the guide field in the ±?/-directions, caused 
by the Lorentz force, and is composed of electrons or 
positrons only. As a result, the guide field creates a sep- 
aration of charges within each island across the layer, 
which in return induces a strong electric field Ey of or- 
der Bq across the layer (see Fig. 13, right panels). We 
also note the induction of an E^ electric field of order 
O.lBo, reversing across the layer, due to the motion of 
the incoming plasma towards the layer at a velocity of 
order /3rccC in the ±y-directions. We also find that the 
strong energy- dependent anisotropy of the particles and 
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Fig. 13. — Spatial distribution of the plasma density (left panels), and electric field strength Ey (right panels) for three different values 
of the guide field strength, B^/Bq = (top), 0.25 (middle), and 1 (bottom), at toji = 221. The figure shows the bottom layer only. The 
plasma density is normalized to the initial drifting particle density inside the layer, no, and the electric field is normalized to Bq. White 
solid lines are magnetic field lines in the xy-plane. 



radiation is preserved with a guide field, although the 
structure of the beam is more complex (high-energy par- 
ticles occasionally point towards high latitudes > 0°). 

The deflection of the particles away from the layer does 
not limit significantly the efficiency of particle accelera- 
tion through the Speiser mechanism. The particle energy 
distribution extends up to 7 sa 2 x 10^. The analysis of 
a sample of particles tracked throughout the simulation 
shows that the energetic particles are inevitably deflected 
away from the layer, but only after several crossings of 
the layer. Fig. 14 shows three typical high-energy parti- 
cles trajectories, projected in the yz-plane. The particles 
are progressively carried outside the layer by the guide 
field, and have time to cool over a Larmor timescale. The 
radiation reaction force gradually overcomes the electric 
force because the pitch angle of the particle to the mag- 
netic field line (and hence the perpendicular magnetic 
field) increases progressively as it exits the layer. In con- 
trast, the ejection of the energetic particles into strong 
B± in the weak guide field case is abrupt: the radia- 
tion reaction force jumps rapidly over a short period of 
time and exceeds the electric force, leading to a rapid 
(sub-Larmor) cooling and the emission of > 160 MeV 
synchrotron photons (Fig. 8). Hence, the presence of 
a strong guide field (S^ > B^) effectively inhibits the 
emission of > 160 MeV synchrotron radiation. Fig. 15 
shows the fraction of the total radiative flux emitted 



above 160 MeV, as a function of the guide field strength, 
from = to = Bq. For B^ = Bq, the > 160 MeV 
fiux is about 50 times smaller than the zero-guide field 
case. We note that the > 160 MeV fiux varies by a factor 
2-3 within our statistical sample of 10 simulations with 
no guide field introduced in Section 4.4, as illustrated in 
Fig. 15. 

A moderate guide field may be beneficial for par- 
ticle acceleration in three-dimensional pair plasma re- 
connection. It tends to suppress the development 
of the drift kink instability that occurs only in 3D, 
which may quench non-thermal particle acceleration 
(Zenitani & Hoshino 2008, see also the discussion in 
Liu et al. 2011; Sironi & Spitkovsky 2011; Kagan et al. 
2012). A large 3D simulation such as those presented 
here in 2D is currently beyond our computational re- 
sources. We leave this issue to a future study. 

5. SOLUTION TO CRAB GAMMA-RAY FLARES? 

In this section, we discuss the implications of our find- 
ings in the context of the Crab gamma-ray flares. Below, 
we review the observational features of the flares, and 
show how reconnection can naturally explain them: 

• Synchrotron radiation > 160 MeV: All flares 
show synchrotron radiation well above the 160 MeV 
burnoff limit, up to 375 MeV during the April 
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Fig. 14. — Three typical high-energy particle trajectories, pro- 
jected onto the yz-pla,nc, in the presence of a Bz = Bo uniform 
guide field. The initial {tuii = 0) and final [tuji = 662) positions 
of the particles arc marked by the green triangles and the red dia- 
monds respectively. The gray band depicts the initial upper-layer 
thickness 25. 
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Fig. 15. — Fraction of the total isotropic radiative flux emitted 
above 160 MeV, as a function of the guide field strength. The fiux 
is time-averaged over the most active phase of reconnection, i.e. 
221 < tuJi < 442. The blue squares represent the result from each 
simulation, where Bz/Bq = 0, 0.1, 0.25, 0.5, and 1. This figure 
also shows the variations of the flux above 160 MeV within the 
statistical sample of 10 identical simulations (see Section 4.4), for 
Bz = 0. 

2011 super-flare (Buehler et al. 2012). We showed 
in Section 4.3 that particles can be accelerated 
well above the classical radiation reaction limit, 
via the relativistic Speiser mechanism deep inside 
the layer, and radiate synchrotron radiation up 
to > 400 MeV. The most energetic particles go 
through a substantial fraction of the electric poten- 
tial drop available in the simulation, i.e. £max ~ 
eA-ccBoLx « 3 PeV, with p,,, = EjBo « 0.3, 
Bo = 5 mG, and = 2.7 light-days. 

• Hard spectrum: The April 2011 flare unambigu- 
ously shows an extra synchrotron component on 
top of the quiescent emission. The spectrum aver- 
aged over the flare is hard, such that dF,^/dv cx 
v~^-^ (Buehler et al. 2012), i.e., consistent with 



an emitting population of pairs distributed with 
a power-law of index p w 1.4. The energy of the 
distribution is then dominated by the highest en- 
ergy particles. We find that a hard spectrum (al- 
though not as hard as observed) is naturally ex- 
pected in our simulations during the brightest pe- 
riods of high-energy emission (periods 8-9-10 in 
Fig. 10), when the photon index above 100 MeV 
is around 1.4, 1.5. A hard spectrum could explain 
why there is no detectable counterpart at lower en- 
ergies (X-rays, optical, IR, radio). 

• Ultra-rapid time variability: The first short 
intra-fiare variability was detected in the 4- 
day-long September 2010 fiare, over a 12-hours 
timescale (Balbo et al. 2011). During the 9- 
day-long April 2011 flare, there are significant 
variations of the flux over periods < 8 hours 
(Buehler et al. 2012). Hence, there are even 
smaller structures in the flaring region emitting 
each spike of high-energy radiation. Along the lines 
of Cerutti et al. (2012b), we flnd in Section 4.4 that 
super-fast time variability of the observed high- 
energy flux is expected, due to the strong inhomo- 
geneity and anisotropy of the most energetic par- 
ticles responsible for the > 100 MeV emission. We 
see a flare on Earth only when the beam of high- 
energy radiation crosses our line of sight. The sym- 
metry of the sub-flare time profile is also consistent 
with observations. We attribute the overall dura- 
tion of the Crab flares (from a few days to a cou- 
ple of weeks) to the rcconncction timescale, i.e., of 
order the light crossing time of the reconnecting 
region, 6 days in our simulations (Fig. 10). We 
associate the intra-flare time variability (hours to 
days) to the light crossing time of magnetic islands 
generated by the tearing instability, < 6 hours in 
the simulations (see also Giannios 2012 for a simi- 
lar interpretation in the context of super-fast TeV 
flares in blazars). Our study also reveals that the 
power-density spectrum of the > 100 MeV light 
curve is expected to be almost flat (white noise. 
Fig. 11). This is consistent with observations that 
show a flattening of the PDS at high frequencies 
(Buehler et al. 2012). 

• Flux/energy correlation: One of the most 
remarkable features of the April 2011 flare is 
the discovery of a clear correlation between 
the > 100 MeV flux and the cut-off en- 
ergy, such that vF^{> 100 McV)oc e+^t 
(Buehler et al. 2012). This result is interpreted by 
Buehler et al. (2012) as the signature of the rapid 
variations of relativistic Dopplcr-boosted emission 
(Lind & Blandford 1985). We discovered in Sec- 
tion 4.4 that reconnection can reproduce such 
a correlation as well, with a power-law index 
corr « -|-3, and mimic the effect of Doppler beam- 
ing as expected in ultra-relativistic reconnection 
(Lyubarsky 2005; Giannios et al. 2009; Giannios 
2012). 

6. SUMMARY 

The main result of this paper is the first discovery 
of particle acceleration above the classical synchrotron 
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radiation reaction limit, and the associated emission of 
> 160 MeV synchrotron radiation, in numerical sim- 
ulations of collisionless pair plasma reconnection. For 
this purpose, we developed a new, parallel, relativistic 
PIC code, called Zeltron, that includes the effects of 
the radiation reaction force on particle dynamics self- 
consistently. We confirm in every detail the expectations 
of earlier analytical works (Kirk 2004; Uzdensky et al. 
2011; Cerutti et al. 2012a): the most energetic particles 
are almost linearly accelerated and confined deep inside 
the reconnection layer, where radiative losses are small 
while the accelerating electric field is strong. If the par- 
ticle stays long enough within the layer, it can be ac- 
celerated above the radiation reaction limit defined in 
the upstream plasma. Then, once the particle eventually 
escapes the layer, it radiates > 160 MeV synchrotron 
radiation within a fraction of a Larmor gyration. The 
acceleration mechanism that emerges from this study is 
remarkably simple and robust. 

Although this work addresses a fundamental question 
in particle acceleration, it is mainly motivated by the 
mystery of the Crab Nebula gamma-ray flares. We have 
shown in this study that all the puzzling aspects of the 
flares are consistent with a reconnection event in the neb- 
ula. In addition to the > 160 MeV synchrotron radia- 
tion, our simulations can explain the ultra-rapid time 
variability of the observed high-energy radiation. The 
emission originates mostly from extremely anisotropic 
and compact bunches of energetic particles. An exter- 
nal observer sees a flare when the beam of energetic 



particles crosses the line of sight. We found that the 
power-spectrum of the light curve is well-fit by a power- 
law, which hardens with the energy band of the radia- 
tion. In addition, we discovered that there is a strong 
positive correlation between the emitted radiative flux 
and the high-energy spectral cut-off, mimicking the ef- 
fect of a relativistic Doppler boost. This correlation is 
also consistent with the observations of the Crab flares. 
A strong guide field, i.e. > Bq, deflects particles out of 
the reconnection layer and tends to suppress the emission 
> 160 MeV synchrotron radiation. Our results support 
the magnetic reconnection scenario at the origin of the 
Crab flares. 
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